rm(list = ls())

library(tidyverse)
library(stargazer)
library(lfe)
library(broom) #for tidy
library(lmtest) # For coeftest()
library(sandwich) # For sandwich()
library(estimatr)
library(texreg)
library(psych)
library(default)
library(sjmisc) #for row_sums
library(quantreg)

options(tibble.print_min = 10)
source("functions.R")

load("./data_output/persons_years_sample.Rdata")


# Quantile regressions ----------------------------------------------------

my_vars = c("cog", "noncog", "consci", "extro", "educ_y")

quantileModels = function(my_data, my_model, quantiles = c(0.2, 0.5, 0.8)){
  # Plot function
  nest_mincer = my_data %>% 
    group_by(vuosi) %>% nest %>% 
    mutate(quant_mincer = map(data, function(y){
      map_dfr(quantiles, ~rq(as.formula(my_model), tau = .x, data = y) %>% 
                tidy(., conf.int = T))}))
  
  nest_mincer_unnest = nest_mincer %>% 
    select(quant_mincer) %>% 
    unnest(cols = c(quant_mincer)) %>% ungroup %>% 
    filter(term %in% my_vars)
  return(nest_mincer_unnest)
}

quantilePlot = function(my_models, my_coef, my_title){
  my_models %>% 
    #filter(synvuosi != 1962) %>%
    filter(term == my_coef) %>% 
    ggplot(aes(x = vuosi, y = estimate, color = as.factor(tau), shape = as.factor(tau), linetype = as.factor(tau))) +
    geom_point(size = 3) +
    geom_line(size = 0.6) +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = as.factor(tau)), alpha = 0.1, color = NA)  +
    theme_bw() +
    scale_color_discrete("Quantile", labels = c("20%", "Median", "80%")) +
    scale_shape_discrete("Quantile", labels = c("20%", "Median", "80%")) +
    scale_linetype_discrete("Quantile", labels = c("20%", "Median", "80%")) +
    labs(x = "Year",
         y = "Estimate",
         title = my_title) +
    scale_y_continuous(breaks = pretty) +
    scale_x_continuous(breaks = 2000:2015, labels = paste0("'", str_sub(2000:2015, 3, 4))) +
    guides(fill = FALSE)
}

quantile_table = quantileModels(full_sample, "inc_prop ~ cog + extro + consci + as.factor(synvuosi)")
quantilePlot(quantile_table, "extro", "Extraversion")
ggsave("/figures/panel_quantreg_extro_cog+extro+consci.pdf", height = 4, width = 7)
quantilePlot(quantile_table, "consci", "Conscientiousness")
ggsave("/figures/panel_quantreg_consci_cog+extro+consci.pdf", height = 4, width = 7)
quantilePlot(quantile_table, "cog", "Cognitive")
ggsave("/figures/panel_quantreg_cog_cog+extro+consci.pdf", height = 4, width = 7)

quantile_table_log = quantileModels(full_sample %>% filter(inc_labor!=0), "inc_log ~ cog + extro + consci + as.factor(synvuosi)")
quantilePlot(quantile_table_log, "extro", "Extraversion")
ggsave("/figures/panel_quantreg_log_extro_cog+extro+consci.pdf", height = 4, width = 7)
quantilePlot(quantile_table_log, "consci", "Conscientiousness")
ggsave("/figures/panel_quantreg_log_consci_cog+extro+consci.pdf", height = 4, width = 7)
quantilePlot(quantile_table_log, "cog", "Cognitive")
ggsave("/figures/panel_quantreg_log_cog_cog+extro+consci.pdf", height = 4, width = 7)

